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ABSTRACT 

RX J1856.5— 3754 is one of the brightest nearby isolated neutron stars, and consid- 
erable observational resources have been devoted to it. However, current models are 
unable to satisfactorily explain the data. We show that our latest models of a thin, 
magnetic, partially ionized hydrogen atmosphere on top of a condensed surface can 
fit the entire spectrum, from X-rays to optical, of RX J1856.5— 3754, within the un- 
certainties. In our simplest model, the best-fit parameters are an interstellar column 
density Nh ~ 1 x 10 20 cm -2 and an emitting area with R°° » 17 km (assuming a 
distance to RX J1856. 5-3754 of 140 pc), temperature T°° » 4.3 x 10 5 K, gravita- 
tional redshift z g ~ 0.22, atmospheric hydrogen column j/h ~ 1 g cm -2 , and magnetic 
field B w (3 — 4) x 10 12 G; the values for the temperature and magnetic field indi- 
cate an effective average over the surface. We also calculate a more realistic model, 
which accounts for magnetic field and temperature variations over the neutron star 
surface as well as general relativistic effects, to determine pulsations; we find there 
exist viewing geometries that produce pulsations near the currently observed limits. 
The origin of the thin atmospheres required to fit the data is an important question, 
and we briefly discuss mechanisms for producing these atmospheres. Our model thus 
represents the most self-consistent picture to date for explaining all the observations 
of RX J1856.5— 3754. 

Key words: stars: atmospheres - stars: individual (RX J1856.5— 3754) - stars: neu- 
tron - X-rays: stars 



1 INTRODUCTION 

Seven candidate isolated, cooling neutron stars (INSs) have 
been identified by the ROSAT All-Sky Survey, of which the 
two brightest are RX J1856. 5-3754 and RX J0720.4-3125 
(see Treves et al. 2000; Pavlov, Zavlin, & Sanwal 2002; 
Kaspi, Roberts, & Harding 2006 for a review). These objects 
share the following properties: (1) high X-ray to optical flux 
ratios of log (fx/ /optical) ~ 4—5.5, (2) soft X-ray spectra that 
are well described by blackbodies with kT ~ 50 — 100 eV, 
(3) relatively steady X-ray flux over long timescales, and (4) 
lack of radio pulsations. 
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For the particular INS RX J1856. 5-3754, single tem- 
perature blackbody fits to the X-ray spectra underpredict 
the optical flux by a factor of ~ 6 — 7 (see Fig. [8}. X-ray 
and optical/UV data can best be fit by two-temperature 
blackbody models with fcT£° = 63 eV, emission size 
R% = 5.1(d/140 pc) kmQ kT™ t = 26 eV, and R™ t = 

1 The most recent determination of the distance to 
RX J1856.5— 3754 is 160 pc (Kaplan et al., in prepara- 
tion). However, the uncertainties in this determination are 
still being examined. Therefore, we continue to use the previ- 
ous estimate of 140(±40) pc from Kaplan, van Kerkwijk, & 
Anderson (2002) since the uncertainty in this previous value 
encompasses both the alternative estimate of 120 pc from Walter 
& Lattimer (2002) and the new value. 
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21.2 (d/140 pc) km (Burwitz et al. 2001, 2003; van Kerkwijk 
& Kulkarni 2001a; Braje & Romani 2002; Drake et al. 2002; 
Pons et al. 2002; see also Pavlov et al. 2002; Triimper et 
al. 2004), where T°° = T cff /(l + z 9 ), R°° = R cm (l + z g ), and 
R eTa is the physical size of the emission region. The gravita- 
tional redshift z g is given by (1 + %) = (1 - 2GM/Rc 2 )~ 1/2 y 
where M and R are the mass and radius of the NS, respec- 
tively. However, the lack of X-ray pulsations (down to the 
1.3% level) puts severe constraints on such two-temperature 
models (Drake et al. 2002; Ransom, Gaensler, & Slane 2002; 
Burwitz et al. 2003). It is possible that the magnetic axis 
is aligned with the spin axis or the hot magnetic pole does 
not cross our line of sight (Braje & Romani 2002). Alter- 
natively, RX J1856.5— 3754 may possess a superstrong mag- 
netic field (B > 10 14 G) and has spun down to a period 
> 10 4 s (Mori & Ruderman 2003), though Toropina, Ro- 
manova, & Lovelace (2006) argue that this last case cannot 
explain the Ha nebula found around RX J1856.5— 3754 (van 
Kerkwijk & Kulkarni 2001b). On the other hand, a single 
uniform temperature is possible if the field is not dipolar 
but small-scale (perhaps due to turbulence at the birth of 
the NS; see, e.g., Bonanno, Urpin, & Belvedere 2005, and 
references therein). 

Even though blackbody spectra fit the data, one ex- 
pects NSs to possess atmospheres of either heavy elements 
(due to debris from the progenitor) or light elements (due 
to gravitational settling or accretion); we note that a mag- 
netized hydrogen atmosphere may provide a consistent ex- 
planation for the broad spectral feature seen in the atmo- 
sphere of RX J0720.4-3125 (Haberl et al. 2004; Kaplan & 
van Kerkwijk 2005). The lack of any significant spectral 
features in the X-ray spectrum argues against a heavy el- 
ement atmosphere (Burwitz et al. 2001, 2003), whereas sin- 
gle temperature hydrogen atmosphere fits overpredict the 
optical flux by a factor of ~ 100 (Pavlov et al. 1996; Pons 
et al. 2002; Burwitz et al. 2003). However, these hydrogen 
atmosphere results are derived using non-magnetic atmo- 
sphere models. Only a few magnetic (fully ionized) hydro- 
gen or iron atmospheres have been considered (e.g., Burwitz 
et al. 2001, 2003), and even these models are not adequate. 
Since kT ~ tens of eV for RX J1856.5-3754 and the ion- 
ization energy of hydrogen at B — 10 12 G is 160 eV, the 
presence of neutral atoms must be accounted for in the mag- 
netic hydrogen atmosphere models; the opacities are suffi- 
ciently different from the fully ionized opacities that they 
can change the atmosphere structure and continuum flux 
(Ho et al. 2003; Potekhin et al. 2004), which can affect fit- 
ting of the observed spectra. 

Another complication in fitting the observational data 
of RX J1856.5-3754 (and RX J0720.4-3125) with hy- 
drogen atmosphere models is that the model spectra are 
harder at high X-ray energies. On the other hand, obser- 
vations of RX J0720.4— 3125 suggest it possesses a dipole 
magnetic field B « 2 x 10 13 G (Kaplan & van Kerkwijk 
2005). It is probable then that RX J0720.4-3125 (and 
possibly RX J1856.5— 3754) is strongly magnetized with 
B ~ 10 13 — 10 14 G, and its high energy emission is soft- 
ened by the effect of vacuum polarization, which can show 
steeper high energy tails (Ho & Lai 2003). Rather than re- 
sorting to a superstrong magnetic field, an alternate possi- 
bility is that there exists a "suppression" of the high energy 
emission. One such mechanism is examined in Motch, Za- 



vlin, & Haberl (2003; see also Zane, Turolla, & Drake 2004), 
specifically, a geometrically thin hydrogen atmosphere at the 
surface that is optically thick to low energy photons and op- 
tically thin to high energy photons. The high energy photons 
that emerge then bear the signature of the lower tempera- 
ture (compared to atmospheres that are optically thick at 
all energies) at the inner boundary layer (usually taken to 
be a blackbody) of the atmosphere model; this leads to a 
softer high energy tail. Motch et al. (2003) find a good fit 
to RX J0720.4— 3125 in this case by using a non-magnetic 
atmosphere model with kT c g = 57 eV, a hydrogen column 
density t/h = 0.16 g cm -2 , and distance of 204 pc, and as- 
suming a M = 1.4M , R = 10 km NS. 

We examine this last possibility by fitting the en- 
tire spectrum of RX J1856.5— 3754 with the latest par- 
tially ionized hydrogen atmosphere models (constructed us- 
ing the opacity and equation of state tables from Potekhin 
& Chabrier 2003) and condensed matter in strong magnetic 
fields (see van Adelsberg et al. 2005) . The goal of the paper is 
to provide a self-consistent picture of RX J1856.5— 3754 that 
resolves the major observational and theoretical inconsisten- 
cies: (1) blackbodies fit the spectrum much better than re- 
alistic atmosphere models, (2) strong upper limits on X-ray 
pulsations suggest RX J1856.5— 3754 may have a largely uni- 
form temperature (and hence magnetic field) over the entire 
NS surface, (3) the inferred emission size from blackbody fits 
are either much smaller or much larger than the canonical 
NS radius of 10 — 12 km. Because of observational uncertain- 
ties (see Section [3} and the computationally tedious task of 
constructing a complete grid of models, we do not attempt 
to prove the uniqueness of our results; rather we try only to 
reproduce the overall spectral energy distribution and argue 
for the plausibility of our model. 

An outline of the paper is as follows. In Section we 
describe the atmosphere model (a thin, magnetized atmo- 
sphere of partially ionized hydrogen with a condensed sur- 
face) and possible ways of producing this atmosphere. A 
brief description of the observations is given in Section [3] 
The model fits to the observational data are shown in Sec- 
tion 2] Possible mechanisms for the creation of thin hydro- 
gen atmospheres and pulsations implied by our model are 
given in Section [S] We summarize and discuss our results in 
Section |S] 



2 MODELS OF NEUTRON STAR 
ATMOSPHERES 

Thermal radiation from a NS is mediated by the atmosphere 
(with scaleheight ~ 1 cm) . In the presence of magnetic fields 
typical of INSs (B > 10 12 G), radiation propagates in two 
polarization modes (see, e.g., Meszaros 1992). Therefore, to 
determine the emission properties of a magnetic atmosphere, 
the radiative transfer equations for the two coupled photon 
polarization modes are solved. The self-consistency of the 
atmosphere model is determined by requiring the fractional 
temperature corrections AT(t)/T(t) <C 1% at each Thom- 
son depth r, deviations from radiative equilibrium <C 1%, 
and deviations from constant flux < 1% (see Ho & Lai 2001; 
Ho et al. 2003; Potekhin et al. 2004, and references therein 
for details on the construction of the atmosphere models). 
We note that the atmosphere models formally have a depen- 
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dence, through hydrostatic balance, on the surface gravity 
g [— (l + z g )GM/R 2 ] and thus the NS radius R; however, the 
emergent spectra do not vary much using different values of 
g around 2 x 10 14 cm s -2 (Pavlov et al. 1995). Nevertheless, 
we construct models using a surface gravity that is consis- 
tent with the inferred radius obtained from the spectral fits 
in SectionS] (g = 1.1 x 10 14 cm s~ 2 with M = 1.4M Q) 
R = 14 km, and z g = 0.2). Also, though our atmosphere 
models can have a magnetic field at an arbitrary angle Ob 
relative to the surface normal, the models considered in Sec- 
tions [2] and [4] have the magnetic field aligned perpendicular 
to the stellar surface (see Section 15.11 for justification and 
other cases). We describe other elements of our atmosphere 
models below. 

2.1 Partially Ionized Atmospheres 

As discussed in Section Q] previous works that attempted 
to fit the spectra of RX J1856.5-3754 and other INSs with 
magnetic hydrogen atmosphere models assume the hydro- 
gen is fully ionized. The temperature obtained using these 
models (or simple blackbodies) are in the range kT°° ~ 
40 — 110 eV. Contrast this with the atomic hydrogen bind- 
ing energies of 160 eV and 310 eV at B = 10 12 G and 10 13 G, 
respectively. Therefore the atmospheric plasma must be par- 
tially ionized. 

Figure Q] illustrates the spectral differences between 
a fully ionized and a partially ionized hydrogen atmo- 
sphere. The atomic fraction is < 10% throughout the at- 
mosphere, where the atomic fraction is the number of H 
atoms with non-destroyed energy levels divided by the to- 
tal number of protons (see Potekhin & Chabrier 2003; Ho 
et al. 2003). Besides the proton cyclotron line at Xb p = 
1966(B/10 12 G) _1 (l + z g ) A, the other features are due 
to bound-bound and bound-free transitions. In particular, 
these are the s = to s = 1 transition at (redshifted) wave- 
length A = 170 A, the s — to s — 2 transition at A = 110 A, 
and the bound-free transition at A = 61 A. The quantum 
number s measures the B-projection of the relative proton- 
to-electron angular momentum (see, e.g., Lai 2001). For a 
moving atom, this projection is not an integral of motion, 
but nonetheless the quantum number s (or m = — s) re- 
mains unambiguous and convenient for numbering discrete 
states of the atom (see Potekhin 1994) . Because of magnetic 
broadening, the features resemble dips rather than ordinary 
spectral lines (see Ho et al. 2003). 

2.2 Thin Atmospheres 

Conventional NS atmosphere models assume the atmosphere 
is geometrically thick enough such that it is optically thick 
at all photon energies (t\ S> 1 for all wavelengths A, where 
the optical depth is t\ — J \\ ds, the extinction coefficient 
is XX, an d the distance traveled by the photon ray is ds); 
thus the observed photons are all created within the atmo- 
sphere layer. The input spectrum (usually taken to be a 
blackbody) at the bottom of the atmosphere is not particu- 
larly important in determining the spectrum seen above the 
atmosphere since photons produced at this innermost layer 
undergo many absorptions/emissions. The observed spec- 
trum is determined by the temperature profile and opaci- 
ties of the atmosphere. For example, atmosphere spectra are 




A (A) 

Figure 1. Spectra of hydrogen atmospheres with B = 4 X 10 12 G, 
and T°° = 4.3 X 10 5 K. The solid line is for a partially ionized 
atmosphere, the dashed line is for a fully ionized atmosphere, and 
the dotted line is for a blackbody. All spectra are redshifted by 
z g = 0.22. 

harder than a blackbody (at the same temperature) at high 
energies as a result of the non-grey opacities (see Fig. [1]); 
the opacities decline with energy so that high energy pho- 
tons emerge from deeper, hotter layers in the atmosphere 
than low energy photons. 

On the other hand, consider an atmosphere that is ge- 
ometrically thinner than described above, such that the at- 
mosphere is optically thin at high energies but is still op- 
tically thick at low energies (jx < 1 for A < Athin and 
t\ > 1 for A > A t hi n ). Thus photons with wavelength 
A < Athin pass through the atmosphere without much atten- 
uation (and their contribution to thermal balance is small 
since most of the energy is emitted at A > Athin in the case 
of RX J1856.5— 3754). If the innermost atmosphere layer (at 
temperature Tthin) emits as a blackbody, then the observed 
spectrum at A < Athin will just be a blackbody spectrum at 
temperature T = Tthm. Motch et al. (2003) showed that a 
"thin" atmosphere can yield a softer high energy spectrum 
than a "thick" atmosphere and used a thin atmosphere spec- 
trum to fit the observations of RX J0720.4-3125. 

The method for constructing a thin atmosphere does 
not differ significantly from that used for a thick atmosphere 
(see Ho & Lai 2001 for details on thick atmosphere construc- 
tion). The difference is that, instead of computing the model 
to large Thomson depths (where t\ 3> 1 for all A), we go to 
a relatively shallow depth at density pthin- At this shallow 
depth, we adjust the inner boundary condition (depending 
on the model, either a blackbody or a condensed surface 
spectrum; see Section 12. 3[) to account for reflections. The 
temperature profile is mainly determined by the radiative 
flux (through radiative equilibrium) and mean opacities. For 
instance, the profile is nearly the same in the thin and thick 
atmosphere models shown in Figure [2] because the bulk of 
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Figure 2. Model atmosphere temperature profiles with B = 4 X 
10 12 G and T°° = 4.3 X 10 5 K. The dotted and long-dashed lines 
are the "thick" atmosphere and "thin" atmosphere with yn = 
1.2 g cm -2 , respectively (see text for details). The short-dashed 
horizontal line indicates the effective temperature (T c ff = 5.3 X 
10 5 K) of the atmosphere model. 



the radiative energy comes from wavelengths for which the 
atmosphere is optically thick; this is clearly seen in Fig- 
ure [3] which plots the atmosphere spectra seen by a distant 
observer. Thus, the difference in the spectra at short wave- 
lengths basically does not affect the mean (Rosseland-like) 
opacities and the resulting temperature profile. The only no- 
ticeable deviation occurs near the condensed surface and is 
due to the boundary condition imposed at this depth. For 
comparison, Figure [4] shows a profile with still lower pthin- 
In this case, the boundary between the optically thick and 
thin portions of the spectrum lies at a longer wavelength. 
Since the integrated flux is more sensitive to the spectral 
difference at these longer wavelengths (being closer to the 
peak of the spectrum), we see a more noticeable difference 
in the temperature profiles. A qualitatively similar result 
is obtained for the same atmosphere thickness but with a 
higher effective temperature, since the peak flux is shifted 
to shorter wavelengths. 



2.3 Condensed Iron Versus Blackbody Emission 

In addition to atmosphere models in which the deepest 
layer of the atmosphere is assumed to be a blackbody, 
we construct (more realistic and self-consistent) models in 
which this layer undergoes a transition from a gaseous at- 
mosphere to a condensed surface. A surface composed of 
iron is a likely end-product of NS formation, and Fe con- 
denses at p » 561 AZ~ 3/5 B% S g cm -3 « 2.35 x 10 4 g cm" 3 
and T < 10 5 5 B% 5 K ~ 5.5 x 10 5 K for the case consid- 
ered here (Lai 2001); note that there is several tens of per- 
cent uncertainty in the condensation temperature (Medin & 
Lai, private communication; see also Medin & Lai 2006a,b). 



10 17 



10 16 



10 15 



7 10 14 
£ 

o 

7 10 13 



BO 

§3 i° 12 



10" 



io 10 



10 9 



..•/; 
" // 
/ I 

I i ; 
i ; 

j i 

thin 
T x<l 



-n-rq i r~ 

w 

V 

optically thick 

— > 



B=4xl0 12 G 
T» = 4.3xl0 5 K 
z g =0.22 



\ 



w 



\ 



\ 



\ 



\ 



blackbody 
-thin atmosphere 
■thick atmosphere 



N 



10 



10 z 10 3 

X (A) 



10 4 



Figure 3. Spectra of hydrogen atmospheres with B = 4 X 10 12 G 
and T°° = 4.3 X 10 5 K. The dotted and long-dashed lines are the 
model spectra using the "thick" atmosphere and "thin" atmo- 
sphere with yn = 1.2 g cm -2 , respectively (see text for details). 
The short-dashed line is for a blackbody with the same temper- 
ature. All spectra are rcdshiftcd by z g = 0.22. The vertical line 
separates the wavelength ranges where the atmosphere is optically 
thin (t\ < 1) and optically thick (t\ > 1). 
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Figure 4. Model atmosphere temperature profiles (top) and 
spectra (bottom) with B = 4x 10 12 G and T°° = 4.3 X 10 5 K. 
The long-dashed and dotted lines are the "thin" atmospheres with 
j/H = 1-2 and 0.12 g cm -2 , respectively. The short-dashed hori- 
zontal line indicates the effective temperature (T g = 5.3 X 10 5 K) 
of the atmosphere model. All spectra are redshifted by z g = 0.22. 
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The hydrogen condensation temperature is also uncertain 
(see, e.g., Lai & Salpeter 1997; Potekhin, Chabrier, & 
Shibanov 1999; Lai 2001); the results of Lai (2001) give 
T ~ 2x 10 4 Bi2 65 K, which suggests much lower temperatures 
than is relevant for our case. The condensed matter surface 
possesses different emission properties than a pure black- 
body (Brinkmann 1980; Turolla, Zane, & Drake 2004; van 
Adelsberg et al. 2005; Perez- Azorin, Miralles, & Pons 2005); 
in particular, features can appear at the plasma and proton 
cyclotron frequenciefl 

We use the calculations of van Adelsberg et al. (2005) 
to determine the input spectrum in our radiative transfer 
calculations of the atmosphere. However, at the tempera- 
ture (Tthin ~ 7 x 10° K) of the condensed layer relevant 
to our thin atmosphere models that fit the spectrum of 
RX J1856.5— 3754, the input spectrum (where t\ < 1) is 
effectively unchanged from a blackbody (since the tempera- 
ture profile is nearly identical, with AT t hi n ~ 3%; see Fig. [2}. 
Thus there are only slight differences in the resulting surface 
spectrum. This is illustrated in Figure [5j where we show the 
emission spectrum from a B = 4 x 10 12 G condensed iron 
surface at T = 7 x 10 5 K and p = 2.35 x 10 4 g cm~ 3 and 
compare this to a blackbody at the same temperature. The 
deviation from a blackbody is smaller at low angles (with 
respect to the magnetic field) of photon propagation 9 and 
increases for increasing 6, as illustrated by the two angles 
9 = 15° and 60° (see van Adelsberg et al. 2005). Thus for 
most angles 6, the condensed surface spectra at short wave- 
lengths (where t\ <C 1, so that this surface is visible to an 
observer above the atmosphere) are virtually identical to a 
blackbody. On the other hand, the atmosphere is optically 
thick at longer wavelengths, where the condensed surface 
spectra deviate from a blackbody; thus the condensed sur- 
face and the spectral features are not visible. We see from 
Figure [3] that the harder spectrum at high energies in the 
"thick" atmosphere becomes much softer in the "thin" at- 
mosphere and takes on a blackbody shape. In contrast, there 
is a negligible difference where the atmosphere is optically 
thick. 



3 OBSERVATIONS & ANALYSIS 

We collect publically available optical, UV, and X-ray data 
on RX J1856.5— 3754. These data have been discussed else- 
where so our treatment will be brief. First, we assemble the 
optical (B- and JJ-band) photometry from the Very Large 
Telescope (VLT) from van Kerkwijk & Kulkarni (2001a) and 
the Hubble Space Telescope (HST) WFPC2 F170W, F300W, 
F450W, and F606W photometry (Walter 2001; Pons et 
al. 2002) as analyzed by van Kerkwijk & Kulkarni (2001a). 



2 van Adelsberg et al. (2005) consider an approximation in their 
treatment of the ion contribution to the dielectric tensor which 
leads to a spectral feature at the proton cyclotron frequency. 
However, because of the uncertainty in this approximation, the 
strength of the feature is not well-determined. Nevertheless, our 
results are not at all strongly dependent on this feature (or the 
input spectrum at these low energies) because the optical depth 
of the atmosphere t\ 2> 1 at the proton cyclotron (and plasma) 
frequency (see text for discussion). 
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Figure 5. Condensed iron surface spectrum (solid line for photon 
propagation direction 9 = 15° and dot-dashed line for 9 = 60°), 
with B = 4x 10 12 G, T = 7 x 10° K, and p = 2.35 x 10 4 g cm" 3 , 
compared to a blackbody (dashed line) with the same temper- 
ature. All spectra are redshifted by z g = 0.22. The vertical line 
separates the wavelength ranges where the atmosphere is optically 
thin (t\ < 1) and optically thick (t^ > 1). 

We also take the optical VLT spectrum from van Kerk- 
wijk & Kulkarni (2001a) and a STIS far-UV spectrum^ 
The spectra are entirely consistent with the photometry as 
calibrated by van Kerkwijk & Kulkarni (2001a), although 
given the limited signal-to-noise ratio of the spectra, we rely 
primarily on the photometry in what follows. We then use 
the Extreme Ultraviolet Explorer (EUVE; Haisch, Bowyer, 
& Malina 1993) data as discussed and reduced by Pons et 
al. (2002) . Finally, we take the RGS spectrum from the 57-ks 
XMM-Newton observation and the 505-ks Chandra LETG 
spectrum that are discussed by Burwitz et al. (2003) . 

A source of uncertainty is that, as mentioned by Bur- 
witz et al. (2003), the RGS and LETG data are not entirely 
consistent in terms of flux calibration: while they have very 
similar shapes (and hence implied temperatures and absorp- 
tions) the radii inferred from blackbody fits differ by as much 
as 10% and the overall flux by as much as 20%. Since the 
LETG fits in Burwitz et al. (2003) are more consistent with 
those of the CCD instruments on XMM-Newton (EPIC-pn 
and EPIC-MOS2) and in our opinion the current low-energy 
calibration of LETG is more reliable, we adjust the flux of 
the RGS data upward by 17% to force agreement with the 
Chandra data. We do not know for certain which calibra- 
tion (if either) is entirely accurate, so some care must be 
taken when interpreting the results at the 10%-20% level. 
Fully reliable calibration or even cross-calibration at the low- 

3 datasets: O5G701010-O5G701050, O5G702010-O5G702050, 
O5G703010-O5G703050, O5G704010-O5G704050, O5G705010- 
O5G705050, O5G751010-O5G751050, and O5G752010- 
O5G752050. 
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energy ends of the Chandra and XMM-Newton responses 
(< 0.2 keV) is not currently available (see, e.g., Kargaltsev 
et al. 2005), and the detailed response of EUVE compared 
to those of Chandra and XMM-Newton is also unknown. 
Therefore, for accuracy in doing the EUV/X-ray fits, we 
concentrate on the LETG data, which are consistent and 
have high-quality calibration. 

We follow the HRC-S/LETG analysis threads "Ob- 
tain Grating Spectra from LETG/HRC-S Data'Q, "Creating 
Higher-order Responses for HRC-S/LETG Spectra'fJ "Cre- 
ate Grating RMFs for HRC Observations'!^, and "Compute 
LETG/HRC-S Grating ARFs'Qand use CIAC0 version 3.2.2 
and CALDB version 3.2.2. We extracted the dispersed events 
and generated response files for orders ±1, ±2, and ±3. After 
fitting the LETG data, we compare the fit results qualita- 
tively with the RGS and EUVE data; the general agreement 
is good, but we do not use them quantitatively. 



4 ATMOSPHERE MODEL FITTING 

Because of data reduction and cross-calibration issues (see 
Section [3]) and possible variations in the interstellar absorp- 
tion abundances (standard abundances are assumed), we 
do not feel that a full fit of the data is justified at this 
time. Therefore we do not fit for all of the parameters in a 
proper sense nor do we perform a rigorous search of param- 
eter space. Instead, we fit for a limited subset of parameters 
while varying others manually. This allows us to control the 
fits in detail and reduce the computational burden of prepar- 
ing a continuous distribution (in B, T e ff, and i/h) of models, 
yet still determine whether our model qualitatively fits the 
data. 

For a given magnetic field and atmosphere thickness, 
we generate partially ionized atmosphere models for a range 
of effective temperatures. (Note that, since the continuum 
opacity of the dominant photon polarization mode decreases 
for increasing magnetic fields, the required thickness yii of 
the atmosphere increases for increasing B.) We then perform 
a x 2 fit in CIA0 to the LETG data over the 10-100 A range 
(0.12-1.2 keV) for the absorption column density Nh [us- 
ing the TBabs absorption model from Xspec (Wilms, Allen, 
& McCray 2000), although other models such as phabs give 
similar results], the temperature T°°, the normalization (pa- 
rameterized by -R°°), and the redshift z g , where we interpo- 
late over T°°. We obtain a good fit, and Table [T] lists the 
best-fit parameters and models; the radius is given assuming 
a distance di4o = d/(140 pc) = 1 (see footnote [1]). Figurc[6] 
shows the confidence regions for temperature and radius; 
the covariance between the other fit parameters are not as 
significant. While we find a range of magnetic fields that 
give acceptable fits, changes in the magnetic field outside 
this range (but still within B = 10 12 - 10 13 G) and atmo- 
sphere thickness lead to worse fits. At B > 10 13 G, spectral 
features due to proton cyclotron and bound species appear 



4 See http://asc.harvard.edu/ciao/threads/spectra_letghrcs/ 

5 See http://asc.harvard.edu/ciao/threads/hrcsletg_orders/ 

6 See http://asc.harvard.edu/ciao/thrcads/mkgrmf_hrcs/ 

' Sec http://asc.harvard.edu/ciao/threads/mkgarf_letghrcs/ 
8 |http:/ /cxc. harvard.edu/ciao/ 



Table 1. Fits to the X-ray Data 



Atmosphere Blackbody 





Model Parameters 




B (10 12 ) 


3 


4 




j/H (g cm -2 ) 


1.2 


1.2 






Fit Results^ 




AT H (10 20 cm" 2 ) 


1.18(2) 


1.30(2) 


0.91(1) 


T°°(10 5 K) 


4.34(3) 


4.34(2) 


7.36(1) 




17 


17 


5.0 


z 9 


0.27(2) 


0.22(2) 




X 2 /dof 


0.85/4268 


0.86/4268 


0.86/4269 



a Numbers in parentheses are 68% confidence limits in the 
last digit(s). 

^ The formal fit uncertainty is < 10%. However, since the ra- 
dius determination depends on the distance, we conservatively 
adopt the current ~ 30% distance uncertainty as our radius 
uncertainty. 



18.2 I I | I | I | I | I | r 




16.8 1 ' 1 ' 1 < 1 < 1 < 1 < 1 

4.36 4.38 4.3 4.32 4.34 4.36 4.38 

T" (10 5 K) 

Figure 6. 1, 2, 3<r confidence regions for the fit parameters T°° 
and R°° of the B = 4 X 10 12 G atmosphere model. The 3 X 10 12 G 
model gives very similar results. The dot indicates the best-fit 
point given in Table [T] 

in the observable soft X-ray range (though they are likely to 
be broadened; see Section f5 ■ 1 [I . which are not seen. 

To further evaluate the quality of this fit, we fit the same 
LETG data with a blackbody. The blackbody fit yields pa- 
rameters (see Table [TJ that are very similar to those derived 
by Burwitz et al. (2003; see Section [T}. A comparison of the 
residuals for the blackbody fit and the B — 4 x 10 12 G atmo- 
sphere model fit is plotted in Figure [7] Given the compara- 
ble Xr (~ 1) we achieve from our blackbody and atmosphere 
model fits (along with the low-energy calibration uncertain- 
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10 20 30 40 50 60 70 80 
A (A) 

Figure 7. Fit residuals for the blackbody and B = 4 X 10 12 G 
atmosphere model given in Table [T] 

ties), we are confident that the atmosphere model describes 
the observations just as well as a blackbody. 

Next, we examine the quality of the fit to the opti- 
cal/UV data, van Kerkwijk & Kulkarni (2001a) showed that 
these data are well fit by a Rayleigh- Jeans power law (Fx oc 
A -4 ) with an extinction Ay = 0.12 ± 0.05 mag. In our fit- 
ting, we try using both the optical extinction implied by the 
X-ray absorption [Ay = AT H /(l-79 x 10 21 cm" 2 ) mag; Pre- 
dehl & Schmitt 1995] and fitting freely for Ay, but we find 
that these fits are too unconstrained and that the value of 
Ay is not sensitively determined by the data (indeed, this is 
reflected in the large uncertainties found by van Kerkwijk & 
Kulkarni 2001a). As a result, we fix Ay to 0.12 mag. We also 
assume a single value for the reddening (Ry = 3.1) and use 
the reddening model of Cardelli, Clayton, & Mathis (1989) 
with corrections from O'Donnell (1994). We find that our 
best-fit model to the X-ray data also produces a A -4 power 
law but underpredicts the optical/UV data by a factor of 
15— 2O9<0. Looking in detail at the highest quality optical 
data point (the HST F606W measurement), we find that it 
is 15% above our model spectra. However, the error budget 
is 3% (photometric uncertainty), 5% (Ay uncertainty), 10% 
(uncertainty in the optical model flux) , and 5% (uncertainty 
in the fitted optical flux due to the X-ray normalization), 
and therefore the 15% disagreement can easily be explained 
by known sources of uncertainty. 

Figure [8] shows the observations of RX J1856.5— 3754. 

9 Comparing monochromatic fluxes derived from photometry 
with the models is not sufficiently accurate for a detailed, quanti- 
tative analysis. A more accurate method would involve convolving 
the models with the filter bandpasses and predicting monochro- 
matic count-rates to compare with the data (e.g., van Kerkwijk & 
Kulkarni 2001a; Kaplan et al. 2003). However, given the assump- 
tions about extinction and reddening and the level of accuracy of 
this analysis, the first approach will suffice here. 




10' 10 2 10 3 10 4 



Wavelength (A) 

Figure 8. Spectrum of RX J1856. 5-3754 from optical to X- 
ray wavelengths. The data points are observations taken from 
various sources. Error bars are one-sigma uncertainties. Optical 
spectra are binned for clarity: STIS data into 30 bins at a resolu- 
tion of 12 A and VLT data into 60 bins at 55 A resolution. The 
solid line is the absorbed (and redshifted by z g = 0.22) atmo- 
sphere model spectrum with B = 4 X 10 12 G, yn = 1.2 g cm -2 , 
T°° = 4.3 X 10 5 K, and R°° = 17 km. The dashed line is the unab- 
sorbed atmosphere model spectrum. The dash-dotted line is the 
(absorbed) blackbody fit to the X-ray spectrum with = 5 km. 

We also overlay the blackbody and our B = 4 x 10 12 G 
atmosphere model spectra with the parameters given in 
Table Q] As discussed in Section [1] the data are gener- 
ally featureless, while the models show spectral features; at 
B ps (3 — 4) x 10 G, the features due to bound species lie 
in the extreme UV to very soft X-ray range and are thus 
hidden by interstellar absorption. Overall, we see that our 
atmosphere model spectra are generally consistent with the 
X-ray and optical/UV data, while a blackbody underpre- 
dicts the optical/UV by a factor of 6—7. 



5 DISCUSSION 

Two important uncertainties in our model for 
RX J1856.5— 3754 are discussed here. The first is flux 
variability as a result of surface magnetic field and tem- 
perature distributions. The second is the creation of thin 
hydrogen atmospheres. 

5.1 Surface Magnetic Field and Temperature 
Variations and Constraints on Pulsations 

A major unexplained aspect of RX J1856.5— 3754 is the 
strong observational limit on the lack of pulsations (< 1.3%). 
As discussed by Braje & Romani (2002), the low pulsations 
could be due to an unfavorable viewing geometry and/or 
a very close alignment between the rotation and magnetic 
axes. Assuming blackbody emission, they find this is pos- 
sible in ~ 5% of viewing geometries (with R = 14 km). In 
fact, they note that if RX J1856.5— 3754 is a normal radio 
pulsar, then the lack of a radio detection may imply low 
X-ray pulsations, since the pulsar beam does not cross our 
line of sight. On the other hand, atmosphere emission can be 
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Table 2. Parameters for Neutron Star Surface 



magnetic latitude 


B 


Teff 


e s 


(deg) 


(10 12 G) 


(10 5 K) 


(deg) 


0—10 


6 


7 





10—40 


5 


6 


30 


40—70 


4 


5 


60 


70—90 


3 


4 


90 




or 






70—90 










highly beamed in the presence of a magnetic field (Shibanov 
et al. 1992; Pavlov et al. 1994) and magnetic field varia- 
tions over the NS surface will induce surface temperature 
variations (Greenstein & Hartke 1983). There arises then 
the questions of whether our single temperature and mag- 
netic field (strength and orientation) atmosphere model fit 
is physically correct and whether it produces stronger pulsa- 
tions than is observed. In response to the former, clearly our 
model spectrum represents an average over the surface and 
the single temperature and magnetic field indicate "mean" 
values. 

In order to determine the strength of pulsations, syn- 
thetic spectra from the whole NS surface must be calculated. 
Such synthetic spectra are necessarily model-dependent (see, 
e.g., Zane et al. 2001; Ho & Lai 2004; Zane & Turolla 2005, 
2006), as the magnetic field and temperature distributions 
over the NS surface are unknown. Therefore, we adopt a 
relatively simple model: we assume the surface is symmet- 
ric (in B and T) about the magnetic equator and divide the 
hemisphere into four magnetic latitudinal regions. We gener- 
ate atmosphere models for each region with the parameters 
given in Table [2] (note that the magnetic field distribution is 
roughly dipolar and Ob is the angle between the magnetic 
field and the surface normal) . Using an analogous formalism 
to that described in Pechenick et al. (1983) and Pavlov & 
Zavlin (2000), we calculate phase- resolved spectra and light 
curves from the whole NS surface (we assume M — 1.4M© 
and R = 14 km). The resulting phase-resolved spectra have 
bound and cyclotron features that are broadened due to the 
surface magnetic field and temperature variations but oth- 
erwise are not significantly different from the "local" spec- 
trum we use in Section [4] thus they would produce a quali- 
tatively similar fit. From the light curves, we obtain energy- 
integrated pulse fractions for the wavelength range 10-85 A 
(the approximate range of the X-ray observations). 

Assuming the magnetic axis is orthogonal to the ro- 
tation axis, f P < 0.01,0.04,0.07 for £ = 10°, 20°, 30°, re- 
spectively, where / P = (F max - F min )/(F max + F min ) is the 
pulse fraction and ( is the angle between the rotation axis 
and the direction to the distant observer. Extending the ar- 
guments of Braje & Romani (2002) to realistic atmosphere 
emission, we thus find there exist viewing geometries for 
which the pulse fraction is at the observed limits, despite 
the relatively large magnetic field (and hence strong beam- 
ing) and temperature variations. We also note the follow- 
ing: (1) An Ha nebula is found around RX J1856.5— 3754; 
if interpreted as a bow shock that is powered by the rota- 
tional energy loss of a NS magnetic dipole, then the spin 
period P » 0.5(B/10 12 G) 1/2 s (van Kerkwijk & Kulka- 
rni 2001b), so that P ~ 1 s. (2) The emission can be 



strongly beamed at higher energies, so that the pulse frac- 
tion is energy-dependent; pulsations may be more evident 
at the higher X-ray energies, although the soft spectrum of 
RX J1856.5— 3754 means that broad energy ranges are likely 
still the most sensitive to pulsations. 

5.2 Creation of Thin Hydrogen Atmospheres 

From the atmosphere model fits in Section [4] (see also 
Fig. (3J) the critical wavelength Athin ~ 30 A gives an at- 
mosphere column j/h ~ 1 — 2 g cm~ 2 . A hydrogen at- 
mosphere with this column density has a total mass of 
A/ H « 1 x 10 13 (i?/10 km) 2 (y H /l g cm" 2 ) g or about 10" 21 
of the mass of a 1.4 Mq NS. Since the diffusion timescale is 
extremely short in a high gravity environment (Alcock & II- 
larionov 1980) , the atmosphere contains the bulk of the total 
hydrogen budget of the NS. The origin of such thin H lay- 
ers is a problem which we now address by briefly discussing 
possible mechanisms for generating thin H atmospheres. We 
leave a more detailed study to future work. 

Accretion at low rates may create a thin hydrogen layer 
on the NS surface. A thin hydrogen layer of j/h = 1 g cm~ 2 
has a total mass of ~ 6 x 10 -21 Mq. Hence the time- 
averaged accretion rate over the age of RX J1856.5— 3754 
(~ 5 x 10 5 yr; see Section |SJ is ~ 0.8g s _1 . For the case 
of RX J0720.4-3125 with y H = 0.16 g cm" 2 , such a mass 
layer may result if the time-averaged accretion rate onto the 
NS is about 0.06 g s" 1 over 10 6 yr. Such small amounts of 
accretion indicate an enormous fine-tuning problem. Thus 
this process is unlikely to be the origin of the thin hydrogen 
envelope. 

We next consider the possibility that the thin hydro- 
gen layer is a remnant from the formation of the NS. In 
this case, Chang & Bildsten (2003, 2004) point out that dif- 
fusive nuclear burning (DNB) is the dominant process for 
determining what happens to hydrogen on the surface of a 
NS. Due to the power law dependence of the column life- 
time T co i on yn, DNB leaves a thin layer of hydrogen, whose 
thickness depends on the thermal history of the NS and the 
underlying composition. We find the column lifetime from a 
full numerical solution to be r co i w 10 10 yr for a magnetized 
envelope with B = 10 12 G and r col « 10 s yr for B = 10 13 
G (assuming M = 1.4 M , R = 14 km, T eff = 5 x 10 5 K, 
and yH = 1 g cm -2 ). The shorter column lifetime for higher 
magnetic fields is due to the effect of the magnetic field on 
the electron equation of state, which increases the hydro- 
gen number density in the burning layer (Chang, Arras, & 
Bildsten 2004) . Thus there is insufficient time at the current 
effective temperature to reduce the surface hydrogen to such 
a thin column, assuming that the age of RX J1856.5— 3754 is 
~5x 10 5 yr. However, DNB was much more effective in the 
past when the NS was hotter. In fact, during the early cool- 
ing history of the NS, DNB would have rapidly consumed all 
the hydrogen on the surface (Chang & Bildsten 2004). This 
leads to another fine-tuning problem: if we fix the lifetime 
of an envelope to some value representing the time the NS 
spends at a particular temperature, we find that the result- 
ing column scales like j/h oc T~g° (Chang & Bildsten 2003); 
therefore, the size of the atmosphere is extremely sensitive 
to the early cooling history of the NS. As a result, DNB ap- 
pears to be a less likely mechanism for producing such thin 
atmospheres. 
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Finally, we briefly examine a self-regulating mechanism 
that is driven by magnetospheric currents and may pro- 
duce thin hydrogen layers on the surface. Since the accel- 
erating potential could be as high as 10 TeV (Arons 1984), 
high energy particles would create electromagnetic cascades 
on impact with the surface (Cheng & Ruderman 1977). 
These cascades result in e/7 dissociation of surface mate- 
rial, analogous to proton spallation of CNO elements in ac- 
creting systems (Bildsten, Salpeter, & Wasserman 1992). 
Protons would be one of the products of this dissocia- 
tion and would rise to the surface due to the rapid dif- 
fusion timescale. Since protons cannot further dissociate 
into stable nuclei, a hydrogen layer is built up to a col- 
umn roughly given by the radiation length of the ultrarel- 
ativistic electrons, beyond which hydrogen can no longer 
be produced and the above mechanism terminates. The ra- 
diation length of ultrarelativistic electrons depends on the 
stopping physics. In the zero magnetic field limit, the stop- 
ping physics is dominated by relativistic bremsstrahlung, 
and the cross-section is a ~ qfZVq, where ap = e 2 /he 
is the fine structure constant, Z is the charge number of 
the nuclei, and ro = e 2 /m c c 2 is the classical electron ra- 
dius (Bethe & Heitler 1934; Heitler 1954). This gives a stop- 
ping column j/stop ~ 3000 g cm~ 2 for hydrogen. Bethe & 
Ashkin (1953) and Tsai (1974) give a more accurate value, 
2/stop ~ 60 g cm -2 , for the stopping column of atomic 
hydrogen; thus it appears that the resulting columns are 
too thick. However, for sufficiently strong magnetic fields, 
the stopping physics of ultrarelativistic electrons is signifi- 
cantly modified, and the dominant stopping mechanism is 
via magneto-Coulomb interactions (Kotov & Kel'ner 1985; 
see also Kotov, Kel'ner, & Bogovalov 1986). In magneto- 
Coulomb stopping, relativistic electrons traveling along field 
lines are kicked up to excited Landau levels via collisions and 
de-excite, radiating photons with energy "/Huje, where TvjJb 
is the Landau cyclotron energy and 7 is the Lorentz factor 
of the incoming electron. Compared to zero-field relativis- 
tic bremsstrahlung, the magneto-Coulomb radiation length 
is smaller by a factor of ~ na-p (Kotov & Kel'ner 1985). 
Applying the correction to the Bethe & Ashkin (1953) re- 
sult (ystop ~ 60 g cm -2 ), we find this gives roughly the re- 
quired thin atmosphere column, j/ s to P ~ 1 g cm' 2 . Though 
extremely suggestive, this mechanism requires a more com- 
plete study and is the subject of future work (see, e.g., 
Thompson & Beloborodov 2005; Beloborodov & Thomp- 
son 2006). 



6 SUMMARY AND CONCLUSIONS 

We have gathered together observations of the isolated neu- 
tron star RX J1856.5— 3754 and compared them to our lat- 
est magnetic, partially ionized hydrogen atmosphere models. 
Prior works showed that the observations were well-fit by 
blackbody spectra. Here we obtain good fits to the overall 
multiwavelength spectrum of RX J1856.5— 3754 using the 
more realistic atmosphere model. In particular, we do not 
overpredict the optical flux obtained by previous works and 
require only a single temperature atmosphere. [Note that 
this single temperature (and magnetic field) serves as an 
average value for the entire surface.] In addition to the neu- 
tron star orientation and viewing geometry, the single tem- 



perature helps to explain the non-detection of pulsations 
thus far. At high X-ray energies, where the atmosphere is 
optically thin, the model spectrum has a "blackbody-like" 
shape due to the emission spectrum of a magnetized, con- 
densed surface beneath the atmosphere. The atmosphere is 
optically thick at lower energies; thus features in the emis- 
sion spectrum of the condensed surface are not visible when 
viewed from above the atmosphere. The "thinness" of the 
atmosphere helps to produce the featureless, blackbody-like 
spectrum seen in the observations. Using a simple prescrip- 
tion for the temperature and magnetic field distributions 
over the neutron star surface, we obtain pulsations at the 
currently observed limits. 

Based on a possible origin within the Upper Scorpius 
OB association, the age of RX J1856.5— 3754 is estimated 
to be about 5 x 10 5 yr (Walter 2001; Walter & Lattimer 2002; 
Kaplan et al. 2002). Our surface temperature determination 
{kT°° = 37 eV) is only a factor of 1.7 below the blackbody 
temperature {kT°° — 63 eV) obtained by previous works 
(see Sections [1] and |4} and therefore does not place much 
stronger constraints on theories of neutron star cooling (see, 
e.g., Page et al. 2004; Yakovlev & Pethick 2004). It may also 
be noteworthy that RX J1856.5— 3754 is one of the cooler 
isolated neutron stars and possibly possesses the lowest mag- 
netic field [B ta (3 — 4) X 10 12 G]; the lower magnetic field 
implies a more uniform surface temperature distribution and 
weaker radiation beaming. 

We show that the production of thin hydrogen layers 
is difficult via accretion or diffusive nuclear burning. Accre- 
tion requires a time-averaged rate that is extremely small 
and fine-tuned. Diffusive nuclear burning has no effect at 
the current effective temperature but is very efficient when 
the star was hotter, such that the column thickness depends 
very sensitively on the early cooling history of the neutron 
star. Hence these scenarios cannot easily explain the thin 
atmosphere columns. On the other hand, it may be pos- 
sible for magnetospheric currents to produce the hydrogen 
atmospheres if the stopping physics of relativistic electrons 
is dominated by magneto-Coulomb interactions. We note 
that younger neutron stars may possess atmospheres com- 
posed of helium rather than hydrogen; however, the opacities 
of bound states of helium in strong magnetic fields is still 
largely unknown (see, e.g., Al-Hujaj & Schmelcher 2003a,b; 
Bezchastnov, Pavlov, & Ventura 1998; Pavlov & Bezchast- 
nov 2005). 

Finally, the emission radius we derive from our atmo- 
sphere model fits is R°° « 17 (d/140 pc) km (although recall 
the distance and flux uncertainties discussed in footnote [1] 
and Section O respectively). Accounting for gravitational 
redshift (z g ~ 0.22), this yields R cm « 14 km. This is much 
larger than the inferred radius obtained by just fitting the X- 
ray data with a blackbody (Rx ~ 5 km; see Sections [1] and 
|3| . As a result, there does not appear to be a need to resort to 
more exotic explanations such as quark or strange stars (e.g., 
Drake et al. 2002; Xu 2003; see Walter 2004; Weber 2005 for 
a review), at least for the case of RX J1856.5— 3754. On 
the other hand, our radius is small compared to the radius 
derived from fitting the optical/UV data (R^, t ^ 21 km). 
For a 1.4 Mq neutron star, the latter implies a low redshift 
(z g ~ 0.12) and very large intrinsic radius (i?o™t ~ 19 km); 
this is ruled out by neutron star equations of state, while our 
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radius i? cm ~ 14 km only requires a standard, stiff equation 
of state (see, e.g., Lattimer & Prakash 2001). 
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